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0. Introduction. 

This project has been motivated by the study of turbulent fluid boundary layers in 
the wall region (Aubry et al. [1])- Numerical investigations of models for the dynamics of 
fluctuations in the boundary layer reveal the presence of intermittent solutions (“bursts”) 
that are persistent over a range of parameter values in the model and that correspond 
to heteroclinic cycles in the model equations. Armbruster et al. [2] concluded, using the 
dynamical systems analysis, that these intermittent solutions of the model are an essential 
feature and not accidental. 

In earlier papers [7, 8, 10, 11, 16] Doedel and Friedman have developed an accurate, 
robust, and systematic method for computing branches of homoclinic and heteToclinic 
orbits. These are orbits of an infinite period connecting two fixed points of an associated 
system of autonomous ordinary differential equations. Homoclinic orbits have been shown 
to play a fundamental role in phenomena such as bursting in biology, chaotic vibrations of 
structures, chaotic oscillations in chemical reactions, etc. Heteroclinic orbits are equally 
important in the understanding of the global behavior of dynamical systems and also in 
the study of wave phenomena in nonlinear parabolic partial differential equations. 

The original goal of this project was to accurately compute heteroclinic cycles and 
to investigate numerically how these cycles evolve as the problem parameters vary in 
a model 4-dimensional system studied in [1] and [2]. This model system is a singular 
perturbation problem. Since liomoclinic and heteroclinic orbits often arise in the context 
of singular perturbation problems, this project has evolved into the development of general 
efficient algorithms for such situations. We developed our algorithms using model 2- and 
3- dimensional problems studied by Deng [12]. 

The organization of the paper is as follows. We first formulate the problem and describe 
some algorithms for computation of branches of homoclinics and heteroclinics. Then we 
describe application of these algorithms to several problems of interest. 


1. Formulation of the Problem and Review of Some Earlier Results. 

In earlier papers [7, 8, 10] Doedel and Friedman have considered the problem of finding 
a branch of solutions of the system of autonomous ordinary differential equations 

<0 A) = 0, «(•), /(•» •) G R", A € R"\ 

' ‘ b ) lim u(f) = uo, lim u(t) = Uj. 

• ►— oo I— *o© 

The method utilizes the linear approximation of the unstable (for t < T-, — TL > 0, large) 
and stable (for t > T+, T+ > 0, large) manifolds, under the assumption that solutions 
of (1.1) decay exponentially to their limits at ±oo. Since every translation of a solution 
of Eq. (1.1) is also a solution, to remove this “phase” indeterminacy we need to add a 
constraint. The equation 

/ OO C-N 

(/(«(<), A) - f{q(t), A 0 )) • -/(u(f), A )dt = 0 

seems to be, computationally, the most appropriate way to do this. It is obtained by 
requiring that the current solution u(t) be as “ close” as possible to the previously computed 
solution q(t) (see [7] for the discussion). Our principle result, Theorem 2 in [10] can be 
summarized as follows: 

Let (q, A 0 ) be a solution of (1.1), (1-2). Assume thain\ = 2 ~(n u + — n ) > 0, where 

and are dimensions of the unstable and stable manifolds of uq and u\ respectively. 
Under appropriate assumptions on f and appropriate transversality conditions , in a neigh- 
borhood of (q, A 0 ) there exists an unique solution branch (u(«s), A(s)), (u(0),A(0)) = (^, A 0 ), 
of (1.1), (1-2) and for sufficiently large — T_, 7+ an unique branch (u(s), A 7 -(s)) of ap- 
proximate solutions. Here s is the continuation (such as pseudo -arclength, employed by 
AUTO) parameter . 

Moreover , we have an error estimate 

(i.3) \\X(s) - At(-)|| r .* + ||u(s) -u r (s)||^( R ) < C^e 2T -'‘° 

for some fiQ > 0 > p\. 

Similar results were obtained in Beyn [3, 4], 

In Friedman and Doedel [11] and Friedman [16] the numerical method and its analysis 
have been extended to the case of center manifolds and higher order approximation of the 
unstable and stable manifolds. 

Continuation algorithm 1. ( [10]). The algorithm is based on the following equations: 
(1.4) ti'(f) - Tf(u{t), A) = 0, 0 < t < 1, 


a) f{u o, A) = 0, 

b) /(«!, A) = 0, 
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( 1 . 6 ) 


a) /u(«o,A)tooi = fioiwo i, rn% 6 R", M0i € R, * - 1, n o, 

b) / U («1, Wli € R”, ^li € R, * = 1) 

a) | io 0l | = 1, z = l)-", n o, 

( 1,T ) b ) |tuii| = 1, i = 1, 


1 ( f(u{t ),\ ) - /(g(0> A °)) ■ /*(«(*). A )/(«(0i A) = 0, 


71 0 710 

а) u(0) = uo + £q 'y ] Coiwpj, cq, = 1, 

i=l 1=1 

(1.9) n, m 

б) u(l) = «i + ei ^ciiioi,-, £4 = 1. 

i=i t'=i 

Here for x,y € R ra x • J/ denotes the inner product in R", and we denote by 
|.| the l 2 norm in R"; we shall keep the same notation for the inner product and the norm 
for finite dimensional Euclidian spaces of dimension other then n. The above equations 
constitute a system of ordinary differential equations subject to constraints. In concrete 
cases the number of constraints (and correspondingly, the number of scalar variables) can 
often be significantly reduced by simple algebra and variations on the normalization equa- 
tions Equation (1.4) is the differential equation. It is obtained by truncating (1.1a) to an 
interval [T_, T+], T_ < 0 < T+, setting T = |T_| + T+ and then scaling the independent 
variable t so that it varies from 0 to 1. The actual period T therefore appears explicitly in 
(1.4), whereas T_ and T + do not appear explicitly in (1.4) - (1.9). There are n\ problem 
parameters, viz., A t , i = l,...,n A . Equation (1.5) defines two fixed points of the vector field. 
In (1.6a) we assume that the Jacobian f u (u o,A) has no distinct real positive eigenvalues 
Mi with eigenvectors w 0l , and n - n 0 real nonpositive eigenvalues. Similarly, in (1.6b) 
we assume that the Jacobian / u (m, A) has m distinct real negative eigenvalues fi u with 
eigenvectors wu, and n-n\ real nonnegative eigenvalues. Under appropriate assumptions 
on /, by the stable manifold theorem the fixed point uq has a (strongly) unstable man- 
ifold of dimensions to which the linear subspace S 0 = Span({ui 0 ,}r=i) 1S tangent at u 0 . 
The fixed point u\ has a (strongly) stable manifold of dimension ni to which the linear 
subspace JJ\ = Span({uii,}"li) is tangent at uj. 

We have w 0 i,wu, «o,«l € R”- Equation (1.9a) then requires that the starting point 
u(0) of the orbit u(t ) lie in the tangent manifold So at “distance” eo from the fixed point 
u 0 . Similarly, equation (1.9b) requires the endpoint u(l) to lie in U\ at “distance” ej from 
the fixed point u\. Finally (1.8) represents the phase condition . 

The unknowns uq and u\ can be eliminated entirely fiom (1-4) (1.9) bj< using 
(1.9). Then (1.4)-(1.8) represent n coupled differential equations subject to n c = 


3 


2n + (n + l)(no + nj) -f 3 constraints, of which (1.9) is an integral condition. In addi- 
tion to the vector function variable u(t) £ R" we have scalar variables 

A £ R” A , eo, £ R, 

(1-10) fi 0 i,c 0i £ R, £ R", no, 

^1i,ci,£R, wu £ R", 2 = 1, nj . 

The total number of scalar variables equals n v = n\+(n + 2)(n 0 + nj)+ 2. Formally we need 
n v = n c — n for a single heteroclinic connection. Usually we are interested in computing an 
entire branch (one dimensional continuum) of orbits, in which case n v = n c — n + 1. This 
requirement is equivalent to setting the number of free problem parameters 

(1-H) n\ = n - (no + ??i) + 2. 


The period T is kept fixed in the continuation. For T large and to and t\ small, each 
solution on the branch represents an approximate heteroclinic connection. If we want to 
increase the period T, then we can replace one of the problem parameters, say Aj, by T 
(see [10] for an example of such a computation). 
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2. New Algorithms for Computation of Homoclinic 
and Heteroclinic Orbits. 

We formulate the algorithms in the situation when the unstable (stable) manifold is 
1-dimensional, while the stable (unstable) manifold can have dimension greater than 1. 
See also Monteiro [13]. In this case the direction along the unstable (stable) manifold is 
locally defined by the eigenvector corresponding to the positive (negative) eigenvalue, while 
it is more difficult to determine what linear combination of eigenvectors defines locally the 
direction along the stable (unstable) manifold. To be specific and without loss of generality, 
we assume that the eigenvector u>£ defines the direction of the one-dimensional unstable 
manifold Wf oc (u 0 ) at the fixed point u 0 . 

Starting orbits can be obtained by using either AUTO itself or some initial value solver. 
In applications we used KAOS [14] and VODE [15]. 

Continuation algorithm 2 (floating boundary algorithm). Eqs. (1-4) (1-9) are 

modified by dropping the equations which define the direction along the stable manifold. 


The equations 

now are: 

(2.1) 

u'(t)-Tf{u(t),\) = 0, 0<t<l, 

(2.2) 

5 

o 

II 

O 

(2.3) 

o,AK = dfutf, R", /*SeR, 

(2.4) 

£ 

©V 

II 

(2.5) 

f (f{u{t),\) - f{q(t),\ Q )) ■ f u {u{t),\)f{u{t),\)dt 

J 0 

(2.6) 

u( 0) = u o + eou»5, 


The unknown u 0 can be eliminated from (2.1)-(2.5) by using (2.6). Then (2.1)-(2.5) 
represent n coupled differential equations subject to n c — 2 n + 2 constraints. In addition 
to the vector function variable u(t) € R" we have scalar variables 

AeR n \ c 0 GR, 

( 2 - 7 ) $ e R, i^gR". 
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The total number of scalar variables equals n v = n\ + n + 2. As in [7] we are interested 
in computing an entire branch (one dimensional continuum) of orbits, in which case the 
number of free problem parameters 


( 2 . 8 ) 


77 A = 1. 


Continuation algorithm 3 (steering vector algorithm). The equations for this 
algorithm are: 

(2.9) 

u'(t) - Tf(u(t ), A) = 0, 0 < t < 1, 

(2.10) 

f{u o, A) = 0, 
/(«i, A) - 0, 

(2.11) 

fv(uo, A)t«o = Mo w o> w o G R", fio € R, 

(2.12) 

\d\ = 1, 
Kl = i, 

(2.13) 

[ (/(«(<)» A ) - /(9(0« A °)) • \)f(u(t),X) dt = 0, 

Jo 

(2.14) 

а) w(0) = u Q + c 0 Wq, 

б) u(l) = + t] d, d £ R n . 


Again the unknowns uq and Ui can be eliminated from (2.9)-(2.13), by using (2.14). 
Then (2.9)— (2. 13), represent n coupled differential equations subject to n c = 3n + 3 
constraints. In addition to the vector function variable u(t) € R" we have n v = n\ + 2n + 3 
scalar variables 

(2-15) A G R n Vu7g,d G R", co,ci,/iS GR- 

Again we are interested in computing an entire branch (one dimensional continuum) of 
orbits, in which case we have 

(2.16) n\ = 1 , and hence n v = 2n + 4. 
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We next give two algorithms for obtaining starting orbits, which can be used in 
conjunction with Continuation algorithms 1-3. 

Starting orbit algorithm 1 (IVP Solver ) 

Step 1 Assume that uq , Wq, and A are given with |iuq| = 1. Initialize the “distance” eo by 
“small” number , such as 0.0001 and compute the initial value 

(2.17) u{Q) = u 0 + c 0 w%. 

This provides an initial value for an initial value solver such as KAOS [ 14 ]- Next solve 
the initial value problem: (2.17), (2.18), 

(2.18) u(t) — Tf(u(t), A) — 0, 0 < t < 1, 
for “large” time T, such as 20 - 100. 

Step 2 Switch to A UTO. Initialize uq, ivq , ft, T and A from step 1. Read the data generated 
by KAOS and interpolate it. This provides an initial orbit u(t) for AUTO . 

Step 3. Perform continuation by one of the Continuation algorithms 1-3 . Note that in 
the case of the Continuation algorithm 1, one first needs to compute the projection of u(l) 
found at Step 1 onto the subspace spanned by w\ t , i = 1,...,77]. 

Starting orbit algorithm 2 (floating boundary). 

Step 1. Initialize the period T by a “small” number , such as 0.01 , and the “ distance ” 
eo by another “small” number, such as 0.0001. Given uq and Wq, initialize the solution 
by a constant: 

(2.19) u (t) = uq + cqivq, 0 < t < 1. 


Step 2 Perform continuation in the direction of increasing T } while all other parameters 
are fixed, using the equations 

(2.20) u{t) — T f(u{t). A) = 0, 0<*<1, 

(2.21) u(0) = Uq + e 0 u?Q. 


Starting orbit algorithm 3 (steering vector ). 

Step 1, Initialize the period T by a “small” number , such as 0.01, and the “distance” 
eo by another “small” number , such as 0.0001. Given uq and Wq, initialize the solution 
by a constant : 


( 2 . 22 ) 


u ( t ) = tio + 0 < t < 1. 


and initialize d , ei from the equations: 


(2.23) 


u\ + t\d = uq + cqWq, 
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Step 2 Perform continuation in the direction of increasing T using the equations 


(2.24) 


u'(t) — Tf(u(t), A) = 0, 0 < t < 1, 


(2.25) 


a) it(0) = n o + eo«'o> 

b ) u(l) = uj -f Ci d, d £ R" 


(2.26) 


\d\ = 1 , 

Kl = i- 


Remark . The integral “phase” condition is removed for the Floating Boundary and the 
Steering Vector algorithms because its purpose is to prevent the sharp peaks from moving. 
In this method the peaks must move to approach the heteroclinic orbit from the constant 
solution. AUTO allows the user to set the accuracy to which the variables and parameters 
are computed in the continuation procedure, by varying the tolerances. In the beginning 
of the continuation one must set high tolerances for the parameters and variables, because 
the constant solution is per se a poor approximation of a heteroclinic orbit. 
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3. A Model 2-d Singular Perturbation Problem. 


3.1 Formulation of the Problem. 

The system of equations given below is a model problem which we have used to develop 
our algorithms in the case of Singular Perturbation Problems, 


x = (2 — z)a(x — 2) + (z + 2)[a(x — .to) + 0(y — yo)] 

( 3 . 1 ) y = (2 - z)[d(b - a)(x - 2)/4 + by] + (c + 2)[-0(x - t 0 ) + a{y - yo)] 
ez = (4 — z 2 ) [z + 2 - m(x -f 2)] — tcz 

From Deng [12], it is known that for the parameter values a = 1,6 = 1.5, c = 2,m = 
1.1845, a = 0.01,0 = 5,x 0 = -0.1, yo = -2,e = 0.01, d = -3.5 the solution is a twisted 
homoclinic orbit, while with the same values of parameters except d = -0.2 the orbit is 
nontwisted. 

Singular perturbation problems are characterized by the appearance of a small param- 
eter such as e = 0.01, which in this case makes the system of ordinary differential equations 
a stiff system. Stiff equations are systems where the magnitude of one eigenvalue of the 
Jacobian is considerably greater than the magnitude of the other eigenvalues. 

In this system of equations, it is known that a hyperbolic fixed point exists at 
wo = (2,0, —2). The eigenvalues of the Jacobian evaluated at (2, 0,-2) are (/ii,/i2, /*3) ~ 
( — 1847,5,3) with associated eigenvectors (itfj , ). The eigenvalue \i \ — —1847 is 
responsible for the stiffness of the system. The eigenvector w j gives the local direction 
of the stable manifold, while some linear combination of (w^w") defines the direction of 
the unstable manifold. 

Computing the homoclinic orbit for the 3-dimensional problem is formidable for 3 
reasons: (1) the system of equations in Eq. (3.1) is stiff. (2) It is known theoretically 
that the direction of the unstable manifold W^ c (u o) is defined only by the eigenvector iC|. 
However for the numerically approximate problem (which we solve) both w™ and w% will 
define the direction of Wfo c (uo) and the linear combination of w™ and telf is unknown. (3) 
There are several heteroclinic orbits near the homoclinic orbit, and the slightest numerical 
instability will displace the homoclinic orbit to one of these heteroclinic orbits. In view 
of these reasons, we decided to compute a 2-dimensional homoclinic orbit, to gain some 
insight into the intricacies of computing the full 3-dimensional orbit. There was also the 
possibility that we could then compute the 3-d orbit by a homotopy from the 2-d orbit. 

For the 2-d problem the stable and unstable manifolds, Wf oc (uo) and W^ c (ao) are 
each one- dimensional. 

We attempt to compute homoclinic orbits for the two dimensional system of equations, 
obtained from Eq. (3.1) by setting xj = 0 and y = 3.59 x 1(T 3 (since it is known that for 
the 3-d system (x,y,z) = (1.99,3.59 x HT 3 , -1.99) is a fixed point): 


(3.2) 


x = (2 - z)a(x -2 ) + {z + 2)[a(x - r 0 ) + 0{y - j/o)] 
i = ((4 - z 2 )[z + 2 - m{x + 2)] - tcz) ft 
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3.2 Computational Results 

We use the following notation: 


u o = (uoi,uo2) 

w i = (^“1,^12) 

= K,^) 

^1 = (^11 1 <^12) 

«0 

uo,w“,ie£,di,u(0),u(l) G R 2 
e o> c i € R 
A = ( e,d ) G R 2 


Fixed Point 

Eigenvector tangent to the unstable manifold at u 0 

Eigenvector tangent to the stable manifold at uq 

Normalized steering vector components connecting u(l) 
and uq 

Distance between no and u(0) 

Distance between uq and w(l) 


To compute a homoclinic orbit for this system of equations we use the steering vector 
algorithm and after noting that for a homoclinic orbit The homoclinic orbit is obtained 
for Eq. (3.1) by a series of steps. 

Step 1 (a): Obtain an initial homoclinic orbit for 


(3.3) 


u' = Tf(u, A), u = {x,z) 


subject to the boundary conditions: 


(3.4) 

and normalization condition, 


(a) u( 0) = uq + eoiv j* 

(b) t/(l) = Uq + C\(l\ 


(3-5) 1^1 = 1 

for a total of 5 boundary conditions. Note that dj = (d n ,di 2 ) G R 2 . Initialize 
T — 0.01, eo - 10 , u o = (1.99,-1.99), = (l.O, — 5.2 x 10 -4 ) and the tolerances 

€ “’ 6 A = 10 _ “. Begin computation of the orbit along the unstable manifold. Perform 
continuation with respect to T,d\\,d\ 2 ,t\ using the steering vector algorithm for obtaining 
initial orbits. We have now reached Fig. (3.1), where (d — —0.2,6 = 0.01, eo = 
1.0 x 10~ 6 ,T = 1.15 x 10 " 2 ,/xf = -1891), (/<“ = 3.99, wj, = 0.99, ’= -5.0 x lO" 4 ), 

(cj = 1.0 x 10 6 ,dn = -0.99,6/i2 = 5.2 x 10" 4 ), (u 0J = 1.99, t< 02 = -1.99). 
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Step 1 (b): Decrease the tolerances so that e u ,e\ = 10 -8 . Perform continuation with 
respect to T, dn, d) 2 , t\ until d\ has approximately the same direction as the eigenvector 
ru®, which defines the direction of the stable manifold IP/o C ( u o) near fi xe< i point. 
The computation reaches Fig. (3.2) where c\ = 0.078. In Fig. (3.2) ( d = —0.2, e = 
0.01, e 0 = 1.0 x 10~ 6 ,r = 3.89, /t* = -1891), (/if = 3.99, = 0.99, w u n = -5.0 x 10~ 4 ), 
(ci = 0.078, dn = -0.42, d i2 = 0.90), (u 0 i = 1.99, u 0 » = -1.99). 

Step 2 (a): Switch from the steering vector algorithm to the eigenvector algorithm 
on the right boundary (stable manifold). The problem is formulated as follows: 

u' = Tf(u, A), u = (x,z) 

( 3 . 6 ) 

with the boundary conditions 


(a) u(0) — uq + 

( b ) u(l) = uq + t\d\ 


eigenvalue problem conditions 


( 3 . 8 ) 

normalization conditions 


(«) 

(6) fid, = rfdi 


( 3 . 9 ) 


(а) Ml = 1 

( б ) |<(||=1 


and the fixed point conditions 


(3.10) f{ u o, A) = 0 

We now have a total of 12 boundary conditions. At this point in the computation, the 
direction of the steering vector is only an approximation to the eigenvector which defines 
the direction of the stable manifold on the right boundary. Therefore, the tolerances for 
the state variables and the parameters should be set rather high when switching from 
the steering vector approximation to the eigenvector approximation. Experimentally we 
found that we needed e U i e A = 1. Perform a few steps of continuation ( NMX = 5) with 
respect to T, e, ei,dn,di2*u , n> ty i2>/ i ij/*i? lt 0ii u 02 ^ or a H parameters. Note the 

presence of the singular perturbation parameter t. After this step of continuation, the 
components d\\,d \2 will be aligned along the direction of the eigenvector, which defines 
the direction of the stable manifold. We have now reached Fig. (3.3) where ei = 0.0752 
and d\ = (d u = -5.91 x 10“ 3 , d \2 = 0.99998), is now aligned along the eigenvector ttij, but 
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the singular perturbation parameter has changed its value to t = 0.0111645. In Fig. (3.3) 
the parameter values are ( d = -0.2, e = 0.0116, to = 1.0 x lO -6 ,! 1 = 3.85, n\ = —1693), 
{Hi = 3 . 99 , = 0 . 99 , u>“ 2 = -5.89 x 10~ 4 ), (ti = 0.07, d u = -5.9 x 10~ 3 ,di 2 = 0.99), 

(uoi = 1.99, U 02 = —1.99). 

Step 2 (b): Decrease the tolerances to e u , = 10 8 and proceed with continuation as 
in Step 2(a). The singular perturbation parameter e increased in value from 0.01 to about 
1.11, whence the orbit obtained was homoclinic. At intermediate values of t\ the orbit is 
not homoclinic. For example for e\ = 0.371, the orbit is not homoclinic. The orbit obtained 
for e\ = 1.117774 is shown in Fig. (3.4), where it is clear that the orbit is homoclinic. The 
parameter values in Fig. (3.4) are ( d = —0.2, e = 1.12, to = 10 x 1 0 — 6 , T = 8.75, //* = —11), 
{n't = 2.85, = 0.99, Wy 2 = -0.079), (cj = 5.3 x 10“ 13 , </n = -0.59, d 12 = 0.80), 
{u 0 i — 1.20, uq2 = —1.71). 

Step 3: Attempt to decrease e back to the value of 0.01. Freeze the period T from 
Step 2(b) and perform continuation with respect toe, eo, ej , d\ j , d\ 2 , //“, n 1 1 ^01 , ^02 

using the same set of boundary conditions described in Step 2(a). A surprising result was 
noted: (1) e returned to a value of 0.011117. (2) eg had started at a value eo = 10 -6 
in Step 1, but on this return journey, it reached a value of eo = —2.0 x 10“ 12 and 
ej = 1.117 x 10 -9 . This final orbit is shown in Fig. (3.5). The parameter values 
in Fig. (3.5) are {d = -0.2, e = 0.011, e 0 = -2.7 x lO” 10 ,? 1 = 8.75, nt = -1693), 
{n't = 3.99, = 0.99,1^, = -5.9 x 10~ 4 ), (e, = 1.1 x 10- 9 ,dn = -5.9 x 10" 3 , d u = 0.99), 
(uoi - 1 .99, 1/02 - —1.99). 

Step 4 : We now attempt to compute a branch of homoclinic orbits with respect to 
{d, e). Add the integral phase condition and formulate the problem as follows: 

/011 , u 1 = Tf(u, A), u = {x,z) 


with the boundary conditions 


( 3 . 12 ) 


(a) u{ 0) = u 0 + e 0 ?e“ 

(b) ti(l) = «o + t\d\ 


eigenvalue problem conditions 


( 3 . 13 ) 

normalization conditions 


(«) />? = 
(b) fid , = 


( 3 , 14 ) 


(«) Ml = 1 

(t) |</i| = 1 


12 



and the fixed point conditions 


(3.15) /(«o,A) = 0 


(3.16) J (/(«(*), A) -/(«(0, A 0 )) -j t S(u(t),\)dt = a 

— OC 

Now perform continuation with respect to e, d, to, ci, dn, d]?, refj , u»“ 2 , /zf, /Zj, uoi, uo2- 
This is now a two parameter continuation problem with respect to the parameters (e,d). 
The parameter d starts at d = —0.2 ancl continues on until d = —8000 without any 
significant change in the orbit. The stiffness of the problem is not altered by this variation 
in the parameter d. 

3.3 Figures 

Fig. 3.1. (d = -0.2, e = 0.01, e 0 = 1.0 x 10~\T = 1.15 x lO" 2 ,/^ = -1891), 
= 3.99, uft = 0.99, tcf 2 = -5.0 x 10“ 4 ), (ci = 1.0 x 10~ 6 ,d„ = -0.99, d 12 = 5.2 x 10~ 4 ), 
(«oi — 1-99, uq 2 = —1.99). 


Fig. 3.2. (d = -0.2, e = 0.01, e 0 = 1.0 x 10 — 6 , T — 3.89, /<f = -1891), (/if = 3.99, wf, = 
0.99, uij *2 = 5.0 x 10 4 ), (ej = 0.078, du = — 0.42, c /12 = 0.90), (izqi = 1.99, U 02 —1.99). 


Fig. 3.3. (d = -0.2, e = 0.0116, e 0 = 1.0 x 10- 6 ,T = 3.85, p\ = -1693), (/if = 
3.99, wj, = 0.99, < 2 = -5.89 x 10~ 4 ), (c, = 0.07, dj, - -5.9 x 10~ 3 ,d 12 = 0.99), 
(uoi = 1.99, uq 2 = -1.99). 


Fig. 3.4. (d = -0.2, c = 1.12, co = 1.0 x lO ” 6 , T = 
0.99, u>f 2 = -0.079), ( Cl = 5.3 x 10- 13 ,d n = -0.59, d 12 


8.75, //f = -11), (/if = 2.85, tof, = 
= 0.80), (uoi = 1.20, U 02 = -171). 


Fig. 3.5. (d = -0.2, e = 0.011, c 0 = -2.7 x 10- ,0 ,T = 8.75, / 1 ? = -1693), (/cf = 
3.99, u> n = 0.99, wf 2 = 5.9 x 10 4 ), (ej = 1.1 x 10 ^,dn - — 5.9 x 10 3 ,dj 2 = 0.99), 

(uoi = 1.99, uq2 = —1.99). 
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4. A Model 3-d Singular Perturbation Problem 

4.1 Formulation of the Problem 

The system of equations given below is the model problem for applying continuation 
algorithms to Singular Perturbation Problems, 

i = (2 — z)a(x — 2) + (z + 2)[cv(t — .To) + 0{y — j/o)] 

(4.1) y = ( 2 - z)[d(b - a)(x - 2)/4 + by] + [z + 2)[-/?(x - t 0 ) + a{y - y 0 )] 
tz = (4 — z 2 ){z + 2 — m(r + 2)] — ccz 

From Deng [12], it is known that for the parameter values a = 1,6 = 1.5,c = 2,m = 

1.1845,a = 0.01,/? = 5, to = — 0.1, yo = — 2,e = 0.01, d = —3.5 the solution is a twisted 
homoclinic orbit, while with the same values of parameters except d = —0.2 the orbit is 
nontwisted. 

Singular perturbation problems are characterized by the appearance of a small param- 
eter such as f = 0.01, which in this case makes the system of ordinary differential equations 
stiff. Stiff equations are systems where the magnitude of one eigenvalue of the Jacobian is 
considerably greater than the magnitude of the other eigenvalues. 

In this system of equations, it is known [12] that there exists a hyperbolic fixed point 
near uq = (2,0, —2) for (e,d) = (0.01,— 0.2). A more accurate solution with the IMSL 
subroutine DNEQNF yields w 0 = (1.99469,3.59524 x 10“ 3 , -1.997886) for Eq. (4.1). The 
eigenvalues of the Jacobian evaluated at « 0 are (/ij, % ( — 1891,3.99,5.99) with asso- 

ciated eigenvectors w* = (-5.3 x 10 -3 , 5.5 x 10 -3 , 0.99), w“ = (0.99, 5.2 x 10 _2 ,-5.2 x 10“ 4 ), 
w 2 = (5-2 x 10 3 , 0.99, — 2.7 x 10~ 6 ). The eigenvalue //, = —1891 is responsible for the 
stiffness of the system. The eigenvector gives the local direction of the stable manifold 
Wf oc { u o), while some linear combination of (tu“,u;“) defines the direction of the unstable 
manifold W£ c {u 0 ). ' ] 

To compute the non-twisted homoclinic orbit, we used the steering vector algo- * 

rithm. The boundary u(0) was displaced from the fixed point u 0 by e 0 along the eigen- : 

vector la, as outlined in [13]; this corresponded to the unstable manifold. The boundary j 

u(l) was attached to the steering vector d\ . The results of the computation are shown 
in the attached Figures. 

4.2 Computational Results. 

Drawing on insight obtained from the two dimensional problem, we attempted to 
compute homoclinic orbits for the three dimensional system : : 


(4.2) u - Tf(u, A), u = (x,y,z) 

which symbolically represent Eq. (4.1). To compute a homoclinic orbit for this system 
of equations, we will go through a number of steps, with each step having a different set 
of boundary conditions. 
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We use the following notation: 


uo : 

= (uoi » 

U02 ) Uo3 ) 


= K 1 

i uj“ 2 , u>“ 3 ) 


= Ki 

, ri>22, u>23) 


= Ki 


dl : 

= (dm 

dn, dn) 


co 

ci 


Fixed Point 

First eigenvector tangent to the unstable manifold at u o 
Second eigenvector tangent to the unstable manifold at uo 
Eigenvector tangent to the stable manifold at uo 
Normalized steering vector components connecting u(l) and 
uo 

Distance between uo and u(0) 

Distance between uq and u(l) 


Step 1 (a): Initialize T = 0.01, e 0 = I0 -7 , (e,d) = (0.01, -0.2) and set the tolerances 
to e u ,e\ — 10 -2 . We attempt to start computing the homoclinic orbit along the unstable 
manifold W? oc (u Q ) from the boundary u(0). W£ c (u 0 ) is defined by some linear combination 
of the eigenvectors, ,w%). This linear combination is not known a priori , so following 
Deng [12] we start the computation of the orbit along the eigenvector w “ . We obtain an 
initial orbit for Eq. (4.1) formulated as follows: 


(4.3) u' = Tf(u, A), u = (x,y,z ) 

with the 6 boundary conditions, 

(a) u(0) = uo + eo{iCi cos k + u >2 sin «:} 

(4 ' 4) (6) u(l) = u 0 + e\d\ 

and normalization condition, 

(4.5) |di|=l 

Note that d\ = (d\i , d^). Perform continuation with respect to T, ei , dy , dj 2 , d\z 

using the steering vector algorithm for obtaining an initial orbit as described in [13]. 
Typically, we use NMX = 5 with NBC = 7, NINT = 0, NT ST = 25, NCOL = 5. 

Remark: At Step 1(a), we assume that the solution to Eq. (4.3) remains constant 
over the interval T = 0.01. The tolerances are set to the high value e u ,e \ = 10 2 in the 
event that the assumption is not true and the solution does vary substantially 

Once the solution reaches Fig. (4.1), AUTO will have computed an accurate enough 
solution, so that the next step of computation can proceed with much lower tolerances. In 
fact, in the next step, all we do is to lower the tolerances. 

Step 1 (b): At this point AUTO has been able to find an initial orbit, so we decrease 
the tolerances to e„, t\ = 10~ 8 . Perform continuation with respect to T, ei,dn,di2>^13, 
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using the same set of boundary conditions as in Step 1(a) until d\ has approximately the 
same direction as the eigenvector w f , which defines the direction of the stable manifold near 
the fixed point. In Fig. (5.2) T = 4.3, t\ = 0.3, dj = (7.5 x 10 -3 ,— 2.9 x 10 -3 ,0.95) and 
this is a fairly good approximation to the eigenvector wf = (—5.3 x 10 — 3 , 5.5 x 10 — 3 , 0.99) , 
which defines the direction of the stable manifold, VV^ c (uo). At this point we still use 
pseudo-arclength continuation with NBC = 1,NINT = 0,NTST — 25, NCOL = 5. We 
now attempt to switch from the steering vector to the eigenvector approximation at the 
boundary u(l). 

Remark : It is important to watch the continuation process very carefully in Step 
1(b) to see when d \ , is approximately aligned along the eigenvector w j. Fig. (4.3a), 
shows what happens when u(l) is “too close” to the fixed point for this stage of the 
continuation (ej = 0.1), but d\ = (0.04, —0.99, —2.2 x 10 -5 ) is a hopeless approximation to 
w\ = (—5.3 x 10“ 3 , 5.5 x 10 -3 , 0.99) . We used the 7 boundary conditions of Step 1(b), and 
attempted to perform continuation with respect to T, k, d \\ , d\o, ^13 in the hope that the 
steering vector d\ would align along the eigenvector for some value of k . Unfortunately, 

this was not the case. Once y(t) undershoots (falls below the fixed point value on the right 
boundary u( 1)) as shown in Fig. (4.3b), it is impossible to compute the homoclinic orbit, 
no matter how close one is to the fixed point. 

Step 2 (a) : Switch from the steering vector algorithm to the eigenvector algorithm 
on the right boundary (stable manifold), i.e. Solve the following equation: 

u = Tf(u, A), u = (x,y,z) 

with the 6 boundary conditions, 

(а) ?*(0) = uo + eo{w| cos k + w •“ sin /c} 

(б) u(l) = Wo + eidi 

9 eigenvalue problem conditions, 


(4.6) 


(4.8) 


(«) />? = ,<?«»? 
TO />; = 

(c) f°J, = Mi* 


3 normalization conditions, 


(4.9) 


(“) l ! »“l = 1 

W Kl = 1 

(c) 1*1=1 


and the 3 fixed point conditions, 
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(4.10) /(»o,A) 

We now have a total of 21 boundary conditions. Continuation is carried out with 

respect to the variables: T,e , (ci, d\\ , c/ 12 , d\z), (/ f 2’ w 2\ ,u; 22’ w 2z)' 

(moi,« 02,«03), /*!,«■ At this point in the computation, the steering vector has only ap- 
proximately the same direction as the eigenvector which defines the direction of the stable 
manifold on the right boundary. Therefore, the tolerances for the state variables an t e 
parameters should be set rather high when switching from the steering vector approxima- 
tion to the eigenvector approximation. Experimentally we found that we needed e u , t\ — 
Note that Eq. (4.8c) forces a change from the steering vector algorithm to the eigenvector 
algorithm. After this step of continuation <1\ is aligned along the direction of the eigenvec- 
tor inf which defines the direction of the stable manifold W, s oc {u 0 ). We have now reached 
Fig. (4.4), where e x = 0.345 and d n = -4.87 x 10“ 3 , c/ i2 = 5-06 x 10 ,d 13 = 1.0 are now 

aligned along the eigenvector twf, but the singular perturbation parameterjias changed its 
value to e = 0.0092. The other parameter values are d = -0.2, e = 9.2 x 10 , T = 4.3, eo - 

10- 7 e, =0.3, ac = 5.7 xlO- 4 , (/if = -2039, d„ = -4.9 xl0- 3 ,d 12 = 5.1x10 ,dn = 1.0), 

(p* L 4, = 0.99, w\ 2 = 0.05, = -4.9 x 10~ 4 ), (/*? = 5.99, u>Ji = 4.8 x 10 ,u>? 2 - 

0.99, w 2Z = -2-3 x 10~ 6 ), (u 0 i = 1.99, u 0 2 = 3.3 x 10 3 ,«o3 = -1-99). 

Remark: It is also possible to switch from the steering vector approximation to the 
eigenvector approximation via a homotopy . 

Step 3 (a): Decrease the tolerances to e U i e A = 1° _1 ’ set NTST = 100, NCOL = 
5 = _1.0 x 10 -4 and use Natural Parameter Continuation to solve: 

u' = r/(u, A), u = {x,y,z) 

(4.11) 

with the integral phase condition, 


(4.12) 


OO 

I (/(u(().A)-/W<).A°)) j t f(n(t),m 


the 6 boundary conditions, 

(a) ti(0) = tio T cos K T w 2 

(4.13) (fr) u(l) = uo + £\d\ 

9 eigenvalue problem conditions, 


(4.14) 


(«) />“ = 

(b) = fiw l 

(c) fid i = p\d x 
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3 normalization conditions, 


(a) K| = 1 

(b) K| = 1 

(c) \di | = 1 


(4.15) 


and the 3 fixed point conditions, 

(4-16) f(uo, A) 

We now have a total of 21 boundary conditions plus an Integral Phase Condi- 
tion. Continuation is carried out with respect to the variables: e\,T,e, (d u , <i 12 , d ]3 ), 

{^\i w \\i w \2i w \z )) (d 3 ’ ^21 > w 22 » ^23 )’ ( w 0l , u 02, u 03), /<f, «,£()• Note that continuation is 
carried out with t\ as the primary parameter, which in conjunction with the negative 
value of DS and Natural P arameter Continuation , ensures that the length of the steering 
vector continuously decreases, pulling the orbit towards the fixed point. Note also that e 0 
and k are allowed to vary so that the correct linear combination of eigenvectors is found on 
the two dimensional unstable manifold, c ( u o ) to compute an accurate homoclinic orbit. 

In Fig. (4.5), the singular perturbation parameter e = 9.33 x 1 0 — 3 , e 2 = 0.344. 
The parameter values are: d = -0.2, e = 9.3 x 10 — 3 , 7" = 4.3, e 0 = 1.9 x 10 -7 ,e] = 
0-34, k = -5 x 10 -6 , (/<f = -2025, dji = -4.9 x 10 _3 ,di 2 = 5.1 x 10~ 3 ,d 13 = 0.99), 
(/if = 4,ia" = 0.99, tt>i 2 = 0.05, uft = -4.9 x lO” 4 ), (/if = 5.99, = 4.9 x 10- 3 ,tof 2 = 
0.99, u> 23 = —2.4 x 10 6 ), (uqi = 1.99, tio 2 — 3.35 x 10 _3 ,i<o 3 = —1.99). The next step in 
the computation is to decrease t\ still more, while simultaneously attempting to increase 
the accuracy of the computation, which involves decreasing the tolerances e u , t\. 

Remark: It was extremely crucial in this step to set NT ST — 100 to carry out the 
continuation procedure. From the phase space plot of Fig. (4.5), it is clear that there 
are several sharp fronts, and this high value of NT ST = 55 is the only way to accurately 
compute the solution with these sharp fronts. For example, a value of NT ST = 55 allowed 
tolerances of only c u ,e\ = l and this was not sufficiently accurate for continuation. 

Step 3(b): The only change from Step (3a) is that we decrease the tolerances and set 
e u, ca = 10 -8 . The continuation proceeds as before with Natural Parameter Continuation 
with respect to the same parameters as in Step (3a). We stop the computation at the 
point shown in Fig. (4.6), where e 0 = 3.7 x 10 -9 ,ei = 1.0 x 10 _7 ,e = 9.3 x 10 _3 ,r = 5.3. 
The parameter values are: d = — 0.2,e = 9.3 x 10- 3 ,T = 5.3, e 0 = 3.7 x lO -9 ,^ = 
1.1 x 10 -7 , k = -2.7 x 10- 7 , (/tf - -2026, dn = -4.9 x 10“ 3 , d r2 = 5.1 x 10~ 3 ,di 3 = 0.99 
), (//f = 4 ,ti7fj = 0.99, wf 2 = 0.05, < 3 = -4.9 x 10~ 4 ), (/if = 5.99,^ = 4.9 x 10“ 3 ,u;f 2 = 
0.99, = 2.4 x 10 ), (uoi = 1-99, t <02 = 3.5 x 10 3 , uo 3 — — 1.99). 

Remark: The extremely small values of eo and ei indicate that the homoclinic orbit 
has been computed very accurately. In practice, an accuracy of eo,ei = 1.0 x 10 -5 should 
suffice. 

Step 4: We now attempt to compute a branch of homoclinic orbits with respect 
to the parameters (d, e) of Eq. (4.1). Deng [12] notes that the homoclinic orbit is 
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twisted for d = -3.5 and non-twisted for d = -0.2. Consequently, we again use Natural 
Parameter Continuation to decrease d from —0.2 to —3.5. We use the same set of boundary 
conditions as in Step (3a) but perform continuation with respect to d, e, eo, (fi , d\\ , di2, ^ 13)5 
(^,^ 1 ,^ 2 , U'“ 3 ), Ox$,wJ 1 ,“>5 2 ,“>5 3 ), (“ 01 , “ 02 , “03), /*?,«• This is now a two parameter 
continuation problem with respect to the parameters (e,d). The parameter d starts at 
d = — 0.2 and continues on until d — — 3.5 with a significant change in the orbit as shown 
in Fig. (4.7). The parameter values are: d = — 3.5, e = 9.3 x 10 3 , eo = 4.7 x 10 9 ,ei = 
1.1 x 10 -7 , k = -1.4 x 10- 5 , {ft\ = -2024, dn = -5.3 x 10' 3 , c / 12 = 5.5 x 10" 3 ,di 3 = 1.0 
), (ft* = 4 = 0.75, wjfe = 0.66, w u n = -3.7 x 10 “ 4 ), (/4 = 5.9 ,w^ = 4.9 x 10 3 ,w u n = 
0 . 99 ,te 2 3 = -2-4 x 10 -6 ), (u 0 i = 1 . 99,« 0 2 = 1-99 x 10 _ 3 ,u 03 = -1-99 ). The stiffness of 
the problem is not altered by this variation in the parameter d. 

Step 5: Deng [12] notes that a homoclinic orbit exists for e = 0.01. Accordingly, 
we tried to increase e from e = 0.0093 (in Fig. 4.6) to e = 0.01 using the boundary 
conditions of Step (2a). We performed continuation with respect to e,eo, (ei,dn,di 2 ,di 3 ), 
(a^“>2 1 »“'m.“ , 23)’ (“01, “ 02 , “ 03 ), /*!,« using both Pseudo -Arclength 
Continuation and Natural Parameter Continuation. The continuation process did not 
converge. 

We also attempted to use the boundary conditions of Step (3a) and perform con- 
tinuation with respect to c,! 1 , (fi, dn, dj> , dn), (d 1 > w n ’ w \ > > w 'h ) , {d2’> w 2i’ w 22i w 23^ 
(U 01 ,“ 02 ,U 03 ), again using both Natural Parameter Continuation and Pseudo- 

Arclength continuation. Once more, there was no convergence. 

A third attempt to perform continuation, which also met with failure was the following. 
At Step (2a) switch from the steering vector approximation to the eigenvector approxima- 
tion on the boundary u( 1) but hold e constant at t = 0.01 and solve the following problem: 

u = Tf{u, A), t* = 

(4.17) 

the 6 boundary conditions, 

(а) u(0) = uq + cos k + iv *2 sin k } 

(б) it(l) = i/o + c ldi 

3 eigenvalue problem conditions, 

(4.19) fid 1 = n\ di 

and the normalization condition, 


(4.20) |<*i| = l 

We attempted to perform continuation with respect to e,T, (fi,dn,di 2 ,di 3 ), f.i j,/c but 
the continuation process did not converge. 

We therefore conclude that it is difficult, if not impossible to compute the homoclinic 
orbit for e = 0.01. 
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4.3 The Methods Which Failed 


1. The Homotopy From Two Dimensions 

The solution of problems such as the Josephson Junction [11], have been carried out 
using a homotopy from a known solution to a simpler problem. Since we had already 
obtained a homoclinic orbit for the simpler two-dimensional problem defined by the system 
of equations: 

i = (2 - z)a{x - 2) + [z + 2)[a(x - x 0 ) + 0{y - j/o)] 

^ ’ ei = (4 - z 2 ) [z + 2- m(x + 2)] - ccz 

we attempted to find a solution to the system, 

X = (2- z)a(x - 2) + {z + 2)[q(.t - x 0 ) + (3{y - j/o)] 

(4.22) y = 7{(2 - z)[d(b - a)(x - 2)/4 + by) + (z + 2)[-/3{x - x 0 ) + a(y - J/o)]} 
ez = (4 — z 2 )\z + 2 — m(x + 2)] — tcz 

We initialized 7 — 0, when Eq. (4.22) reduces to Eq. (4.21) and having found a homoclinic 
orbit for the system of ODEs in Eq. (4.21), we attempted to increase 7 from 7 = 0 to 
7 = 1, which would give a homoclinic orbit in three dimensions. The advantage of dealing 
with a two-dimensional problem first, is that the stable and unstable manifolds are one- 
dimensional, so we do not have to deal with unknown linear combinations of eigenvectors, 
which is necessary to solve the three dimensional problem. However, this approach did 
not meet with any success; with hindsight one can see that there is a loop in the three- 
dimensional orbit shown in Fig. (4.6), which does not exist in the two-dimensional orbit 
of Fig. (3.5), and this is probably why the homotopy fails. 

2. Increasing the Singular Perturbation Parameter t 

Apart from the fact, that the linear combination of eigenvectors on the two-dimensional 
unstable manifold, Wfc (uq) is unknown, this problem is compounded by the fact that the 
system of ODEs is stiff, which makes the problem of computing homoclinic orbits even 
more difficult. For the two-dimensional system defined by Eq. (4.21), it was found that |/i|| 
could be reduced (thus making the system of ODEs less stiff) by increasing the singular 
perturbation parameter, e. With this “insight”, we attempted to increase t using the 
boundary conditions of Step. (2a) and carrying out continuation with respect to e, 7 1 , 

(ei,dli,</i2,C? 13 ), w 13)’ (/*2 ’ W 21 ’ w 22 > w 2'i ) ’ (“01, «02, "(tt), /*!,«■ 

This method failed and the explanation is as follows: It is known theoretically from 
Deng [12] that a homoclinic orbit exists for e = 0.01, which means that the orbit must 
intersect the stable manifold Wfc (iiq) defined by the eigenvector at the right boundary 
u( 1). For larger values of 6, there is no certainty that the orbit will intersect this 
one dimensional stable manifold; in fact, in general the orbit will not intersect the one 
dimensional stable manifold. To guarantee that the orbit intersects the unstable manifold, 
the manifold would have to be two-dimensional which is not the case for this problem. 
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3. The Initial Value Problem Solver VODE 


VODE is an IVP solver designed specifically to tackle systems of stiff ODEs. We 
intended to compute an initial orbit with VODE for the system in Eq. (4.1), and then 
perform continuation using AUTO, with respect to the parameters (e,d). VODE was 
not able to produce the orbits, which were achieved by AUTO in Fig. (4.6). We used 
the data from Fig. (4.6), namely, e = 9.3 x 10 -3 ,d = —0.2, uoi = 1.995, uo 2 = 
3.55 x 10~ 3 , uq 3 - 1-998, ic“, = 0.99, = 5.21 x 10“ 4 , w\ 3 = -4.92 x 10 -4 , urfr = 

4.92 xlO -3 , u> 2 2 = 0.99, u>“ 3 = -2.4 x 10 -6 , to = 3.05 x 10 -8 , k = -1.89 x 10 -6 . The results 
of the computation with VODE are shown in Fig. (4.8 al-cl), and for comparison the 
results with AUTO are shown in Fig. (4.8 a2— c2). VODE is clearly not able to reproduce 
the sharp fronts, which can be computed with AUTO and worse still, the computation of 
y(t) is hopelessly inaccurate. 

4.4 Figures. 

Fig. 4.1. Initial solution showing x(t) which is constant. y(t) and z(t) are also constant. 


Fig. 4.2. T = 4.3, ei = 0.3, dj = (7.5 x 10" 3 ,-2.9 x 10“ 3 ,0.95) 


Fig. 4.3a. A spurious solution plotted in the phase plane 


Fig. 4.3b. Graph of y(t ) showing undershooting of the y component. 


Fig. 4.4. The other parameter values are d = — 0.2, e = 9.2 x 10 3 ,T = 4.3, eo = 
10" 7 ,ei = 0.3, k — 5.7 x 10 — 4 , (pf = -2039, d u = -4.9 x 10- 3 ,d ]2 = 5.1 x 10“ 3 , c/i 3 = 1.0), 
(/i“ = 4, toft = 0.99, tCj 2 = 0.05,1^ = -4.9 x 10" 4 ), (/x? = 5.99 , = 4.8 x lO" 3 ,^ = 
0.99, rw 2 3 = -2-3 x 10 -6 ), (zxqi = 1.99, tx 02 = 3.3 x 10- 3 ,tx o3 = -1.99). 


Fig. 4.5. The parameter values are: d — — 0.2, e = 9.3 x 10 3 ,T = 4.3, eo = 1-9 x 
10~ 7 ,ei = 0.34, k = — 5x 10~ 6 , (/xf = -2025, d n = -4.9x 10~ 3 , d n = 5.1 x 10" 3 , d n = 0.99), 
(/x“ = 4,^ = 0.99, < 2 = 0.05, «)j 3 = -4.9 x 10“ 4 ), (/x.“ = 5.99, wj, = 4.9 x lO" 3 ,^ = 
0.99, w% 3 = -2.4 x 10" 6 "), (uoi = 1.99,no2 = 3.35 x 10- 3 ,u 03 = -1.99). 
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Fig. 4.6. eo = 3.7 x 10 -9 , ei = 1.0 x 10 7 ,e = 9.3x10 3 ,T = 5.3. The parameter values 
are: d = -0.2, e = 9.3 x 10~\T = 5.3, e 0 = 3.7 x 10~ 9 ,ei = 1.1 x 10“ 7 ,« = -2.7 x 10“ 7 , 
(|if = —2026, d\\ = -4.9 x 10“ 3 , r/i 2 = 5.1 x 10~ 3 , <i 13 = 0.99 ), (/t? = 4,^ = 0.99, u>? 2 = 
0.05, = -4.9 x 10~ 4 ), (/<? = 5.99,11)?! = 4.9 x 10- 3 ,w? 2 = 0.99, = -2.4 x 10" 6 ), 

(uoi = 1.99, uq2 — 3.5 x 10- 3 ,u 03 — —1.99). 


Fig. 4.7. The parameter values are: d = — 3.5, t = 9.3 x 10 3 , eo = 4.7 x 10 9 ,ei = 
1.1 x 10 ~ 7 ,k = —1.4 x 10“ 5 , (/if = — 2024, dn = —5.3 x 10 — 3 , r /] 2 = 5.5 x 10~ 3 , c /13 = TO 
), {n\ = 4, u)jj = 0.75, u>j 2 = 0.66, = -3.7 x 10 " 4 ), = 5.9, = 4.9 x 10 " 3 ,u ;? 2 = 

0.99, u )? 3 = —2.4 x 10 ~ 6 ), (uqi = 1 . 99,^02 = T99 x 10 _ 3 ,uo 3 = —1.99 ). 


Fig. 4.8 Comparison of results obtained with VODE (al-cl) and AUTO (a2— c2). The 
parameter values are t = 9.3 x 10~ 3 ,d = —0.2, uoi = 1.995, U 02 = 3.55 x 10 — 3 , 1103 — 1.998, 
wfr = 0.99, u)? 2 = 5.21 x 10" 4 , iw? 3 = -4.92 x 10~ 4 , w u n - 4.92 x 10“ 3 , iv% 2 = 0.99, te? 3 = 
-2.4 x 10~ 6 , e 0 = 3.05 x lO” 8 , « = -1.89 x 10~ 6 . 
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5. The Fluid Mechanics Problem 


5.1 Parameter Assignment in AUTO 

For the generic four dimensional problem, the parameter assignment is as follows 


«o 

e odi^o(4) 

^ 0 . , WQi 

e 0i 

Ui 

eid,^i(4) 

ei» 


par(101)-par(104) left fixed point 
par(106)-par(110) left steering vector 
par(lll)-par(130) 4 left eigenvalues/ vectors 
par(131)-pa.r(135) 4 left eigenvector weights 
par(151)-par(154) right fixed point 
par(156)-par(160) right steering vector 
par(161)-par(lS0) 4 right eigenvalues/ vectors 
par(181)-par(185) 4 right eigenvector weights 


6 

en 

ei2 

6 

C21 

«22 

•S34 


par(l) 

par(2) 

par(3) 

par(4) 

par(5) 

par(6) 

par(7) sign for Eqs. (3), (4) 
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5.2 Formulation of the Problem 

We study this system of equations which are discussed in Eq. (2.5) in the paper by 
Armbruster et al [2]: 


(a) i\ = X\X 2 + 2/12/2 + ^1 (6 + e\\r] + ei 2 r|) 

( b ) yi = xiy-2 - yiX-2 + yi (6 + e u rj + e n r\) 

(5.1) (c) X 2 = ±(.rj - y\ ) + X- 2 (£2 + e 2 i?’i + e 22 r\) 

( d ) i/2 = ± 2 xiyi + y> (£> -f e 2 \r\ + 622^) 

2 2,2 j 2 2,2 

r 1 =^ 1 + ?/i and r 1 = x 2 + y 2 

This system of equations corresponds directly to Eq.(28) of Aubry [1]: 

V‘2 — C 4j _ 0(^4^ + W4IV2) + v 2{ a 2 + ^ 22 r 2 + ^ 24^4 ) 

W 2 = C 4j _ >(^2^4 — V 4 W 2 ) + w 2 { a 2 + d22 r 2 + ^24 r 4 ) 

(5.2) f; 4 = C2,2(^2 _ w i) + v *{ a 4 + d<Yir\ + d 44 r|) 

W 4 = 2 c 2 ;>V 2 W> + uq(a 4 + (Iy 2 ^\ + ^44 r 4 ) 

2 2,2 J 2 2,2 

r 2 = v 2 + w 2 and r \ — v 4 + w 4 

The systems in Eq. (5.1) and Eq. (5.2) are obtained from the Navier-Stokes equations [1]. 
The fluctuating component of the velocity is expanded as a Fourier series in the spanwise 
and streamwise directions. A Galerkin projection is applied to convert the system of PDEs 
into a system of ODEs. The series is then truncated to retain only the first few terms in 
the Fourier expansion becuase the Galerkin approximation minimizes the error due to 
truncation. The important parameters in the Armbruster system are and £>> which 
correspond to the parameters a? and a 4 of the Aubry system, ao an d a 4 are related to the 
Heisenberg parameters a\ and ao and the Reynolds’ number Re t by the equation: 

( . cik = o l k + (1 + Q J i / Rer)^l 

(5- 3 ) j , 2 

c k\k-k' - c k\k-k' + a 2 c k\k-k’ 

The Heisenberg parameters aq and a> may be adjusted upward and downward to 
simulate greater and smaller energy losses to the unresolved modes, corresponding to the 
presence of a greater or smaller intensity of smaller-scale turbulence in the neighborhood 
of the wall. This might correspond to the environment just before or just after a bursting 
event which produces a large burst of small scale turbulence which is then diffused to the 
outer part of the layer. 

We wish to investigate the dynamical behavior of the system defined by Eq. (5.1). 
The initial value used was directed along u?q with cq — 10~ 7 as follows: a{t) — uq + cqico. 
Continuation was carried out with respect to the period which was initially set at T = 0.01 
and increased till T = 130. 
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The initial values at the left fixed point uq were the following: 


(5.4) 




(° 

.18 

0 

0 

°\ 

^) — 


0 

0 

0.34 

0 

0 

-0.40 

0 

0 



\ 

0 

0 

0 

0 ) 

«>o 

= (1 0 

0 

0) 




Co 

= 0.18 






u 0 

41 

0 

II 

// 2 / c 22 

) 1/2 ,o) 

1 



5.3 Computational Results 

We use the following notation: 


uo = (^01, ^02, ^03, W 04 ) 
m = (1/11,1/12,^13,1^14) 

^0 

Mo 

^1 = (^11^12^13,^14) 


eo 


Fixed Point near left boundary 
Fixed Point near right boundary 
Eigenvector defining unstable manifold at uq 
E igenvalue corresponding to u>q 

Normalized steering vector components connecting 
1/1 and n(l) 

Distance between uy and w(0) 

Distance between u\ and t/(l) 


For the system of equations defined by Eq. (5.1) with the negative sign used in Eq. 
(5.1 c-d), we attempt to generate Fig. (5c) in the paper by Armbruster [2] for which we 
used the parameter values: en = — 4,ej2 = — l,eoi = — l,e*22 — — 2, fi = —0.03, {2 = 0.2. 
The unstable manifold W(^ c (uo) at uo is one dimensional and its direction is defined by 
the eigenvector w o- The stable manifold W* (u\) at u\ is two dimensional and since it is 
difficult to determine the linear combination of eigenvectors which determine its direction 
we use the steering vector d\ at this boundary. 

Step 1. Initialize the period T by a “small” number, such as 0.01, and the “ distance ” 
eo by another “small” member, such as 10 -7 . Given i/y and with |iuq| = 1 initialize 
the solution by a constant: 

(5.5) u (t) = «o + eou>o> 0 < t < 1. 

and set the tolerances e u = e\ = 10~ 8 and set the number of subintervals and col- 
location points to NT ST = 25 and NCOL — 4 respectively. Eq. (5.6) and Eq. 
(5.7) represent a total of 9 boundary conditions. We perform continuation with respect 
to (T, ei,dn,di2,di3,di4). We have now reached Fig. (5.1), where (T = 16$, n = 
7.5 x lO'Vn = 0.04, d \2 = 0,di3 = 0.99, d u = 0 ). 
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Step 2 We will now attempt to repeat the results of Fig. (5e) in Armbruster (2j, 
where (4i = 0.135, £2 — 0.2). For these parameter values, modulated travelling waves 
coexist with the heteroclinic orbit. Moreover it is precisely at these parameter val- 
ues that a bifurcation occurs from the 2-dimensional heteroclinic orbit shown in Fig. 
(5-1) ( where (£1 = —0.03, £2 = 0.2)^ to a full dimensional heteroclinic orbit. To reach 
(£2 = 0.135,42 = 0.2) we follow a multistep continuation procedure, where the problem is 
formulated as follows: 


(5.6) 


u'{t)-Tf{u{t),X) = 0, 0<t < 1, 


(«) /(«o, A) = 0 
(6) /(«!, A)=0, 


a) u(0) = tio + 

b) u( 1 ) = tti + e\d\, d\ G R” 


(5.9) 


fu{u Oi A)'U)q 


U U 

Po w o> 


U 


e R”, f‘u e R, 


(5.10) 


(«)|di| = i 

m ki = 1 . 


(5.11) 



(/(u(O.A) -/(«((), A')) 


•/,(«(*), -!)/(“( i)A)dt = 0. 


There are 22 boundary conditions plus an Integral Condition. Continuation is performed 
with respect to the 20 parameters. 

Step (2a): Perform continuation with respect to (4i ,£>,T), (uoi> « 02 , «03, ^ 04 ), 

(un, U12, U13, U14), (dn , d\2, ^13, d]4 ), (fig, Wqj, u>q2, W03, w 04 )• The continuation process 
fails to converge at the terminal values of (4i = —3.9 x 10 3 , ^2 — 0.198), 

(T = 124, eo = 10 -7 ,ei = 10~ 6 )), («oi = 0,u O 2 = 0,tfo3 = 0.315, u 04 = 0), (u u = 0,ui 2 = 

0 ,ui 3 = —0.315, U 14 — 0), (dn = 0.14, di2 = —1.4 x 10 ^,di 3 = 0 . 98 , di 4 = 0.065 ), 

(Mo = 0.2, < = l.ufo = 0 ,< 3 = 0,wg 4 = 0 ). 

Step (2b): Using these terminal values, perform continuation with respect to 

(4li42, c o), (noi,U02,U03,U04), («11 , «12, « 13, ?i 14), (dll , di2, d]3, d i4 ), (/io , ^oi , ^’o2i ^)3i ^04 )• 

The continuation fails to converge at the terminal values (4i = 0.1,42 = 0.21), 
(T= 124,c 0 = 8.5 x 10 _9 ,ei = 10 -6 )), (u 0 i = 0, « 0 2 = 0,u 0 3 = 0.323, u 0 4 = 0), (uu = 
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0,1112 = 0,ui3 = -0.323, u 14 = 0), {du = 0.79, d n = -6.7 x 10 Via = 0.019, d u = 0.060 
), (//“ = 0.326, ufo = 1,^02 = 0^03 = 6^04 = 0 )• 

Step (2c): Using these terminal values perform continuation with respect to 

(6>6> c l)» («01,«02U‘03,W04), («11,«12,W13,U14), (^11,^12,^13,^14 ), (/*(} i u, 01 » W 02’ ^03’ ^04 

). Continuation proceeds till we reach Fig. (5.2) where (£1 = 0.138,^2 — 0.21), 
(r = 124 ,e 0 = 8.5 x 10 -9 ,ei = 2.3 x 10 -6 )), («oi = 0,iz 02 = 0,u O3 = 0.323, 1x04 = 0), 
(un = 0, 1/12 = 0,ui3 = -0.323 ,um = 0), {d u = 0.96, d l2 = -8.1 x 10 ,<*i3 = 
4.5 x I0~\dn = 0.268 ), (/ig = 0.326, = l,u>g 2 = 0,u>J 3 = 0,u?ff 4 = 0 ). 

Step 3: We will also attempt to reproduce the bifurcation diagram for the branch of 

heteroclinic orbits, shown in Fig. 3 in Armbruster [2]. Using the same set of boundary con- 
ditions as in Step 2, increase the number of subintervals to NTST = 55. Using the initial 
values from Fig. (5.1), perform continuation with respect to (6,6, T), (“01, “02, “03, “ 04 ), 
(mi, 1X12, 1113 , U14), (^11.^12,^13 , du ), (/«!,<!, <1^04 )■ The bifurcation diagram of 
(6,6) is shown in Fig. (5.3), where the initial values are (6 = -0.03,6 - 0.2) and the 
final values are (6 = -0.036,6 — 0.09). The entire bifucation diagram in Fig. (5.3) has a 
number of bifurcating branches. We have traced out only one branch. This may account 
for the fact that our results do not match those of Armbruster [2]. 

5.4 Figures 

Fig. 5.1 (T = 168, ei = 7.5 x 10 _9 ,c/n = 0.04, d n = 0,</i 3 = 0.99, d u = 0 ). 


Fig. 5.2 (6 = 0.138,6 = 0.21), (T = 124, e 0 = 8.5 x 10 -9 , ei = 2.3 x 10 6 )), («oi = 
0,u 02 = 0,ao3 = 0.323, u 04 = 0), (un = 0,u 12 = 0,u 13 = -0.323, u u = 0), (d u = 
0.96, d i2 = -8.1 x 10“ 7 ,(/ 13 = 4.5 x 10 _3 ,f/ H = 0.268 ), (/1 q = 0.326, uigj = l,ti>g 3 = 
0,< 3 = 0 ,<= 0 ). 

Fig. 5.3 The initial values are (6 = —0.03, 6 — 0.2) and the final values are 

(6 = -0.036,6 = 0.09). 
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6. Conclusions and Recommendations. 

Homoclinic and heteroclinic orbits are orbits of an infinite period connecting two fixed 
points of an associated system of autonomous ordinary differential equations. Homoclinic 
orbits have been shown to play a fundamental role in phenomena such as bursting in 
biology, chaotic vibrations of structures, chaotic oscillations in chemical reactions, etc. 
Heteroclinic orbits are equally important in the understanding of the global behavior of 
dynamical systems, turbulence, and also in the study of wave phenomena in nonlinear 
parabolic partial differential equations. 

In earlier papers Doedel and Friedman have developed an accurate, robust, and 
systematic numerical method and derived error estimates for the computation of branches 
of homoclinic and heteroclinic orbits. The idea of the method is to reduce a boundary 
value problem on the real line to a boundary value problem on a finite interval by using a 
local (linear or higher order) approximation of the stable and unstable manifolds and then 
study the reduced problem using a continuation software package such as AUTO. 

Theoretical analysis of homoclinic and heteroclinic orbits is often conducted in the con- 
text of singular perturbation problems. In this paper we have refined and extended ealier 
algorithms of Doedel and Friedman using 2 model singular perturbation problems and a 
turbulent fluid boundary layers in the wall region problem. We have thus considerably 
extended the range of applicability of our algorithms 
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